function [erdat,hp1,hp2]=er_results(tsdat,stepdat,fitdat,titlestr,core,nsteps,Aeff)
%function [erdat,hp1,hp2]=er_results(tsdat,stepdat,fitdat,titlestr,core,nsteps,Aeff)
%
% modified from erode_results2.m originally by Steve Suttles
% 
% modified by Patrick Dickhudt - 20-Dec-2010
% - removed evals
% - changed so any mass eroded prior to assigned initial critical stress is
% subtracted.

%calc values needed to find tau_c
 disp(['Core ' core ])
    disp(' ')

    figure
    [ax,h1,h2] = plotyy(tsdat.etime,tsdat.tss,tsdat.etime,tsdat.tau);
    set([h1 h2],'linewidth',1.5)
    ylabel(ax(1),'TSS (g l^-^1)')
    ylabel(ax(2),'Applied Stress (Pa)')
    xlabel(ax(1),'Elapsed time (s)')
    xlim(ax(1),[min(tsdat.etime) max(tsdat.etime)])
    xlim(ax(2),[min(tsdat.etime) max(tsdat.etime)])
    
    title([titlestr ' TSS and applied stress time series'])
    
tau_cin=input('Enter Initial Critical Shear Stress Value, in [Pa] if increase at 0.1 use 0.075: ')

for i1=1:nsteps

    lamb(i1) = fitdat(i1).lambda;
    
    if lamb(i1) <= 1e-9
        lamb(i1)=1e-9;
    end    
    
    tau_app(i1) = mean(tsdat.tau(stepdat.idx.(['step' num2str(i1)])));
end
    
ersteps = find(tau_app >= tau_cin);  % find steps were erosion occured as defined user    

e_ce_i = cumsum(stepdat.em)./Aeff; % cumulative sediment mass collected (kg m-2)

% calculate and subtract mass collected prior to beginning of erosion - may
% be associated with underestimate of background TSS or may represent
% material in suspension overlying core prior to flushing
em_corr = e_ce_i(min(ersteps)-1);  
e_ce = zeros(1,length(e_ce_i));
e_ce(ersteps) = e_ce_i(ersteps) - em_corr;
e_cb = [0 e_ce(1:end-1)];

E_0 = nan(1,nsteps);
E_0_part1 = (e_ce(ersteps)-e_cb(ersteps)).*lamb(ersteps);
E_0_part2 = 1-exp(-lamb(ersteps).*stepdat.tstep(ersteps));
E_0(ersteps) = E_0_part1 ./ E_0_part2;

%start routine to find critical shear stress tau_c0

   
    tau_e = nan(1,max(nsteps));
    tau_ce = tau_cin.*ones(1,max(nsteps));
    
for i2 = ersteps
    tau_e(i2)=tau_app(i2)-tau_ce(i2-1);
    tau_ce(i2)=tau_ce(i2-1)+(1-exp(-lamb(i2).*stepdat.tstep(i2))).*(tau_e(i2));
end

    M=E_0./tau_e;

f1=figure;
clf

hp1(1)=plot(tau_ce,e_ce,'-ob');
set(gca,'ydir','reverse')
set(gca,'XaxisLocation','top')
title(['Core ' titlestr])
xlabel('{\tau}_c , [Pa]')
ylabel('Eroded Mass, [kg ^. m^{-2}]')
ylim=get(gca,'ylim');
xlim([0 0.6])
namestr = ['tau_c' ' core ' titlestr]
set(gcf,'name',namestr)
set(gca,'linewidth',1.5)
set(hp1,'linewidth',1.5)
set(hp1,'MarkerEdgeColor','k')
set(hp1,'MarkerFaceColor','b')


f2=figure

hp2=plot(M,(e_cb+e_ce)./2,'square-r')
set(gca,'ydir','reverse')
set(gca,'XaxisLocation','top')
title(['Core ' titlestr ])
xlabel('M, [kg ^.  m^{-2}^. s^{-1} ^. Pa^{-1}]')
ylabel('Eroded Mass, [kg ^. m^{-2}]')
set(gca,'xticklabel',get(gca,'xtick'))
set(gcf,'name',['Erosion Rate Constant, M,  Core' titlestr ])
set(gca,'linewidth',1.5)
set(hp2,'linewidth',1.5)
set(hp2,'MarkerFaceColor','r')
set(hp2,'MarkerEdgeColor','k')


%save figures

saveas(f1,['tau_' core],'fig')
saveas(f2,['M_' core],'fig')

%make table of results in structured variable and output to csv file

erdat = struct('core',core,'step',stepdat.stepn,'tau_app',tau_app,...
    'e_cb',e_cb,'e_ce',e_ce,'lambda',lamb,'tau_ce',tau_ce,'tau_e',...
    tau_e,'E_0',E_0,'M',M,'Aeff',Aeff)
    
erdat.fields = {'erdat = erosion data processing results';
                'core = core name';
                'step = step number';
                'tau_app = average applied stress (Pa) calculated from RPM time series';
                'e_cb = cumlulative eroded mass at beginning of step (kg m-2) - this is mass collected minus background TSS and any other material collected prior to onset of erosion';
                'e_ce = cumulative eroded mass at end of step (kg m-2) - this is mass collected minus background TSS and any other material collected prior to onset of erosion';
                'lambda = exponential decay constant for each step (s-1), (if < 1E-9, set to 1E-9)';
                'tau_ce = critical stress at end of step (Pa)';
                'tau_e = excess stress at beginning of step, taub - tauc (Pa)';
                'E_0 = erosion rate at beginning of step, calculated from fit to time series (kg m-2 s-1);'
                'M = erosion rate "constant" (kg m-2 s-1 Pa-1)';
                'Aeff = effective surface area of core (m2)'};
              
fid = fopen([core '_results.csv'],'w')

fprintf(fid,'step,tau_applied,time,e_cb,e_ce,lambda,tau_ce,tau_e,E_0,M\n')
fprintf(fid,'%g,%g,%g,%g,%g,%g,%g,%g,%g,%g\n',[erdat.step erdat.tau_app' stepdat.tstep' ...
    erdat.e_cb' erdat.e_ce' erdat.lambda' erdat.tau_ce' erdat.tau_e' erdat.E_0' erdat.M']')
    
fclose(fid)


